Imaging Sky Background Estimation

This Jupyter notebook provides some tips on estimating the sky background for relatively complex scenes. It also demonstrates ways to evaluate the quality of this sky estimation.

Introduction

Sky estimation is one of the tricker aspects of image processing. This is in part because the "sky" is part of the astronomical scene. Some of what is considered background may be in the foreground (thermal background in the detector, scattered light, or zodiacal light). Sometimes the objects of interest overlap (galaxies in front of other galaxies). This notebook does not address the "de-blending" problem of overlapping galaxies, but does outline some of the techniques for dealing with large-scale patterns in the sky.

The Photutils manual has an extensive discussion of background estimation, which is the basis for much of what is in this notebook.

Imports

Matplotlib setup for plotting

There are two versions

Create an image with a somewhat pathological background pattern

This image has a pixel-to-pixel RMS of 0.1, and the background non-uniformities we add are of comparable level to the pixel-to-pixel noise. So you can't get rid of the non-uniformities without doing a fair amount of smoothing. We'll add:

This particular test case was deviously chosen to be one where just adding higher orders to (say) a Chebyshev polynomial surface fit is likely to do a poor job. Similarly, increasing the order of a bicubic spline is not likely to be satisfactory.

Set up grid and the random number seed for the simulation

The test is best illustrated with a 2000x2000 image, but shrink_factor can be used to make the notebook run fast for testing. The random number seed is set to allow for repeatability.

Create the nasty sky background

Blemishes are inserted as sinc-functions with random centers and widths and amplitudes.

Add a gradient on top of the blemishes

Look at it with noise added and then smoothed a bit

Smooth with a 10x10 boxcar kernel

Add some random sources

We'll create a catalog of sources that have elliptical Gaussian profiles with a variety of fluxes, sizes, and position angles. The sizes are kept relatively small relative to the structure in the background. Store these sources in an Astropy table for convenience.

Add these sources to the image to create a scene.

Routine for making a plot to compare two images side by side

Try using Photutils Background2D mesh grid as a first attempt

Now we need to make a first-cut estimate of the background. Photutils provides a flexible interface to estimating the background in a mesh of cells. We use the Background2D routine here, with a fairly typical set of parameters. See the documentation for more information.

This tasks creates a mesh of rectangular cells covering the image. The unmasked pixels within each cell are "averaged" using a configurable robust statistic. These averages can then be "averaged" on some larger scale using a configurable robust statistic. This filtered set of averages is then used to feed an interpolation routine to make a smooth background. The default interpolation is a bicubic spline, but we will illustrate inverse-distance weighting interpolation later on in the notebook.

Try adjusting the arguments to Background2D to see the effect.

Make a couple masks. The make_source_mask routine has a few options. Try changing them to see what they do. The aim here is to find and mask sources without seeing an appreciable enhancement of sources in regions where the background is high. There are several parameters that can be adjusted:

Try reducing the mesh size of the background estimation in the Background2D task. The smaller mesh size definitely helps, but there is a hint in the residual image on the right that sources are biasing the background estimate.

Ring median filtering the image

It's clear we are going to do more iteration to do the source removal. Before doing that, let's try another approach to removing the sources: the ring-median filter. To do this, create a filter kernel that is larger than the sources of interest. Use the scipy median-filtering routine to slide this across the image, taking the median of all the pixels within the ring. This basically gives a local-background estimate for each pixel. (While we don't illustrate it here, it's possible to use scipy.ndimage.generic_filter and pass it a function to compute something other than a median.)

The cell below sets up the kernel.

Compare the scene (minus the mean background) to the filtered scene.

Subtract the ring-median filtered image from the scene as the first cut at background subtraction. Convolve this with a kernel to smooth for source detection, and threshold that to identify pixels that are part of sources.

It's often useful to grow the mask. This can be accomplished by convolving with a filter. Here, we adopt a circular tophat filter.

Mask the image and see how many sources are still remaining

Let's make a fancier plot that shows background residuals, the residuals with the sources and source mask overlayed, and the background-subtracted scene.

Try this fancier plot. The default scaling of the first two plots has a harder stretch in the colormap, so we can now see that there are ring-shaped residuals around all the bright sources. These are not readily apparent in the third plot, which has the stretch we used for the earlier plots.

Set up some routines for convenience in iterating the mask

Next we will try iterating in fitting the background and progressively removing sources at lower and lower S/N. We'll want to inspect at each step. Here are some functions to reduce typing:

Estimate the background on a finer scale after masking sources

Try a different interpolation algorithm.

The Shephard inverse-distance weighting searches for the N grid points nearest to the pixel of interest and weights them as $1/(D^p + \lambda)$ where $D$ is the distance to the neighbor, $p$ is a power, and $\lambda$ is a parameter to make the weighting of the neighbors smoother closer to the pixel of interest.

Make a deeper source mask.

This runs three iterations, with different kernel widths and different tophat sizes for growing the mask. Try varying sigma, filter_fwhm and tophat_size to see how they affect the sources. The aim is to mask sources that are easily visible, but leave plenty of pixels for tracing the background.

Redo the background estimation using this new mask.

Play with n_neighbors, box_size, filter_size and exclude_percentile to see how they affect the residuals.

Check for bias in due to the galaxies

Since we are dealing with a simulation, where we know truth, we can evaluate the biases directly. (In the case of a real image, we can't do that; we'll suggest another test further down).

Tabulate the residual background under the central 3x3 pixels under each galaxy vs the flux of the galaxy, separately for those that are masked and those that are unmasked. First we need to get these values. Take the mean for 3x3 pixels centered on each source. Keep track of the ones that were masked in fitting the background and the ones that weren't so we can check separately for any bias.

Make a list of random positions and do the same measurement, as a control.

Make lists of the masked and unmasked sources

Define a function fit_line to fit a line to the data points, and another function fit_and_plot to plot the data points and the best-fit line together.

Plot the residuals vs. source flux, separately for the masked, unmasked, and randomized source positions. Make the sizes of the markers proportional to the fluxes of the sources. Ideally, this should be centered at 0.0 for all three cases, with no correlation with source flux. Of course, it is nearly impossible not to have a bias where the sources were not masked, since they then contribute to the background estimate. In this particular case, there is a small offset in the residuals when the source was not masked, and evidence of a trend with the flux of the source.

Plot the residuals vs. source radius, separately for the masked, unmasked, and randomized source positions. The sizes of the markers are proportional to the fluxes of the sources.

Testing for bias on real data

With real data, we can't take the residual relative to noiseless sky. However, we can check for evidence that the background under the masked areas is statistically higher than the background in the random areas. Given that astronomical sources have wings and these can't be subtracted without modeling them, it is very hard to achieve perfection here. This test will tell you the magnitude of the potential bias (in flux per pixel) and whether or not it looks significant.

Routines to facilitate looking at the RMS of the residual background as a function of scale

One way to evaluate whether the sky-subtraction looks sensible is to see whether the RMS is dropping sensibly as a function of scale. It should drop linearly with the area size of the image sampled.

To do this test, we would like to look at unmasked regions, but have them all have the same S/N. So we need a way to set up a mesh grid that has no masked pixels in each of the mesh areas.

Calculate the statistics for the different background estimates. These are the residuals of the masked "scene" -- i.e. with the sources in it, but masked as well as one might do for a real scene. These are to be compared to the "perfect" case of the noisy sky background minus the noiseless sky background, where the statistics are computed in the same unmasked cells.

Plot the results relative to the perfect case (we've commented out the initial bkg1 estimate just to make the scale more visible for the better estimates). It is worth noting that the bkg1 is typical of what one might get from SExtractor, which offers only a single pass for the background estimation.

Space Telescope Logo